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A quantum phase transition from paramagnetic to ferro- 
magnetic phase is driven by a time-dependent external mag- 
netic field. For any rate of the transition the evolution is 
non-adiabatic and finite density of defects is excited in the 
ferromagnetic state. The density of excitations has only log- 
arithmic dependence on the transition rate. This is much 
weaker than any usual power law scaling predicted for pure 
systems by the Kibble-Zurek mechanism. 

PACS numbers: 03.65.-w, 73.43.Nq, 03.75.Lm, 
32.80.Bx, 05.70.Fh 

According to the book [1] our understanding of quan- 
tum phase transitions is based on two prototypical mod- 
els: the quantum Ising chain and the Bose-Hubbard 
model. It is well established in the exactly solvable 
Ising model [5-7] and there are indications in the Bose- 
Hubbard model [8] that essential properties of a dynami- 
cal quantum phase transition are captured by the Kibble- 
Zurek (KZ) mechanism of defect formation [2]. This the- 
ory was originally developed for thermodynamic transi- 
tions and tested both by numerical simulations [3] and by 
experiments in finite-temperature condensed matter sys- 
tems [4]. It is based on the observation that, because of 
the critical slowing down, a system driven across its crit- 
ical point must go out of thermal equilibrium no matter 
how slow is the transition rate. Using only combination 
of general casuality and universality arguments it shows 
that in a transition from a disordered to ordered phase 
the system ends in a non-equilibrium state with finite 
ordered domains of average size 
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Here z and v are critical exponents and tq is time of 
the transition (or quench). No matter how slow is the 
transition the system does not have enough time to order 
throughout its whole volume and average size of ordered 
domains is limited to a finite £ which is a power of the 
transition time tq. 

The short list of the two prototypical models must 
be supplemented by a prototypical disordered system 
i.e. the random Ising chain defined by the Hamiltonian 
[10,11] 



N 



71 = 1 



with periodic boundary conditions 
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Here J ra 's are random ferromagnetic couplings, J n > 0, 
and h is external magnetic field. This model has a quan- 
tum critical point at h c = exp(lnJ„) separating ferro- 
magnetic (h < h c ) from paramagnetic (h > h c ) phase. 
When h — > oo the ground state becomes fully polarized 
along the x-axis, 

I— —...), (4) 

but when h = there are two degenerate ferromagnetic 
ground states 
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The critical point at h c is surrounded by the Griffith 
regime (or phase) of infinite linear susceptibility. 

The random Ising model is important because real con- 
densed matter spin systems are random. Their random- 
ness is not a mere small perturbation on top of pure 
models because the disorder is changing their universality 
class. For example, the randomness of J n in the model 
(2) is changing its universality class with respect to the 
pure Ising chain with a constant J n = J. No matter 
how weak is the randomness of J n renormalization group 
transformations drive the model towards an infinite dis- 
order fixed point [11]. As a result, the random Ising 
chain has v = 2 instead of v = 1 and, more importantly, 
z — > oo instead of z = 1. The diverging dynamical ex- 
ponent marks a qualitative difference between the pure 
and the disordered model. As the dynamical exponent 
is relevant to dynamical phase transitions, its singularity 
is suggesting that an outcome of the transition in a dis- 
ordered system can be qualitatively different from that 
in its pure counterpart. This expectation has never been 
tested in the theory of dynamical phase transitions. 

In the model (2) a dynamical paramagnet-ferromagnet 
transition is not adiabatic and, as a result, the state after 
transition is a superposition over excited states like 
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with finite ferromagnetic domains separated by domain 
walls or kinks. The diverging dynamical exponent z im- 
plies that, in first approximation, the size of the ordered 
domains £ in Eq.(l) does not depend on the quench time. 
No matter how slow is the transition the average density 
of kinks d ~ £ remains the same. 

In second approximation we can follow general princi- 
ples of the KZ mechanism [2] to derive £ when z — > oo. 
These ideas were originally proposed for classical tran- 
sitions [2], but they were generalized recently to quan- 
tum phase transitions [8,9,5,6]. The basic idea is that 
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when a system is driven through a quantum phase tran- 
sition, then the evolution of its state is adiabatic when 
the system is away from the critical point, but it must 
be non-adiabatic in a neighborhood of the critical point 
where the gap A between the ground state and the first 
excited state tends to zero. It is convenient to measure 
the distance from the critical point by a dimcnsionlcss 
parameter e = hc ^ h ■ The evolution is non-adiabatic 
between — e and +e, where ±e are the two points when 
transition rate equals the gap A. On one hand, assuming 
that e(t) can be linearized near e = 0, 



e(t) ~ — , 
Q 



we can estimate the rate as 
hand, we know that for |e| 

a ~ kr 



(7) 



On the other 



the gap behaves like 

, (8) 



see Ref. [11]. The rate and the gap equal at e when 
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Here a = 0(1) is a single fitting parameter. Quick anal- 
ysis of this equation shows that e — + when tq — > oo: for 
very slow transitions the evolution is non-adiabatic only 
in a very close neighborhood of the critical point, where 
we can use Eq.(8) self-consistently. 

In the adiabatic-impulse approximation the state of the 
system follows its instantaneous ground state before — e, 
then its evolution becomes impulse between — e and +e 
when the state does not change because reaction of the 
system is too slow as compared to the transition rate, 
and finally the evolution becomes adiabatic again after 
+e. Accuracy of the adiabatic-impulse approximation 
was tested in the quantum context in Refs. [7,9]. The 
state of the system at — e is the ground state with a finite 
correlation length 
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This state does not change until +e when the evolution 
becomes adiabatic again. In this way the correlation 
length £ becomes imprinted on the state of the system 
after the transition. This length becomes the average 
size of ordered domains in the final ferromagnetic state 
(6). For example, in the random Ising model [y = 2) 
approximate solution of Eq.(9) in the limit of tq — > oo 
together with Eq.(10) gives 



(11) 



when In tq 3> 1. Again, as in our first approximation, 
the dependence on tq is very weak as compared to any 
power law scaling. 



As mentioned before, validity of Eq.(8) is limited only 
to close vicinity of the critical point. Further away from 
this point we can expect spectral properties of the model 
to be more like in the pure Ising model. This must be true 
at least for weak disorder, when the random J„'s have a 
narrow distribution around a finite J. As a result, for 
moderately slow transitions, when e is large enough, we 
can expect 
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like in the pure Ising model with v = z = 1 [5,6], see 
Eq.(l). 

Before proceeding with numerical simulations of a dy- 
namical transition we need better understanding of spec- 
tral properties of the random Ising chain. To this end 
I assume for convenience that N is even and make the 
Jordan- Wigner transformation, 
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introducing fermionic operators c„, which satisfy anti- 
commutation relations {c m ,ct} = S mn and {c m ,c„} = 
{c^jcjj} = 0. The Hamiltonian (2) becomes 



H = P+ H+ P+ + P- H~ P- 



where 



P ± = 



N 



i±n< 
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1 ± J] (1 - 24 c 
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(15) 
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are projectors on subspaces with even (+) and odd (— ) 
numbers of c-quasiparticles and 



± N ( 

n=l ^ 
h.C. 
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are corresponding reduced Hamiltonians. The c„'s in H ~ 
satisfy periodic boundary conditions cn+i = ci, but the 
c„'s in H + must obey cn+i = — ci, what I call "antiperi- 
odic" boundary conditions. 

The parity of the number of c-quasiparticlcs is a good 
quantum number and the ground state has even parity 
for any value of h. Assuming that the quench begins in 
the ground state we can confine to the subspace of even 
parity. In this subspace the quadratic H + is diagonalized 
by a Bogoliubov transformation 
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m—1 
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The index m numbers Bogoliubov modes which are eigen- 
modes with positive u> of the stationary Bogoliubov-de 
Gcnncs equations 
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LOU n = 

2hu n - j n {u n+ i + v n+ i) + j n -i {v n -i - u n -i) , 

U)V n = 

-2hu n + j n (V n+ i +U n+ i) + j n -i (V n -i - U n -l) . (19) 

Here I suppressed the mode number m and denned 
jn = Jn for n < N and j^ = — Jn to take into ac- 
count the antiperiodic boundary conditions. The eigen- 
states (u nm ,v nm ) with positive energy oj m > 0, normal- 
ized SO that J2n (\ u nm\ 2 + \v nm I ) — 1 ) define quasipar- 
ticle operators j m = u* nm c n + v n mC n - Each positive en- 
ergy eigenstate has a partner {u~ m ,v~ m ) = {v nm ,u nm ) 
with negative energy — oj m which defines a quasipar- 
ticle operator 7m = {u mn )*c n + v~ m c\ = 7^. Af- 
ter the transformation (18) the Hamiltonian H + = 
(iLlm - ljl m ) equivalent to 

N ( 1\ 

ff+ =E £ ™ hllrn -g) (20) 
m=l ^ ' 

which is a simple-looking sum of quasiparticles. How- 
ever, thanks to the projection P + H + P + in Eq.(15) 
only states with even numbers of quasiparticles belong 
to the spectrum of H. 




Chain size N 

FIG. 1. Average logarithm of the critical gap lnA c as a 
function of N. Each circle is an average over 500 realizations 
of J n and the solid line is the best fit In A c = 0.74-0.58ViV. 

Infinite Ising chain would have zero energy gap at the 
critical point but in numerical simulations we must use 
a finite chain with a finite gap A c at h = h c . As the fi- 
nite gap can alter the adiabaticity condition (9) for slow 
quenches, it is important to know how the gap depends 
on the system size N. For example, in the pure Ising 
chain A c decays only like 1/N and relatively large N 
would be required for accurate numerical simulations. 
In the random Ising chain the "critical gap", which is 
a sum of energies of the two lowest energy quasiparticles 
A c = uji + 0J2 at h = h c , can be found by diagonaliza- 
tion of Eqs.(19). The best fit to an average over different 
realizations of J n is 



lnA c = 0.74 - 0.58VA , (21) 

see Figure 1. For a given tq (logarithm of) the critical 
gap is much less than (logarithm of) the transition rate 
~ Tq 1 provided that N 3> In 2 (tq). This condition can 
be easily met on a relatively small lattice. It is not quite 
surprising that the condition is equivalent to £ <C N, 
compare Eq.(ll). 

Finally, after these preparations, I can proceed with 
numerical simulations of a dynamical quantum phase 
transition. In my simulations I assumed that the inde- 
pendent random J„'s have a uniform distribution in the 
range (0,2) and the critical field is h c = cxpln J n = |. 
For the sake of simplicity the transition is driven by a 
linear quench 

h{t) = - — h c (22) 
T Q 

with t G (—00, 0). The system is initially prepared in the 
ground state at a large initial value of the magnetic field 
hi i.e. in the Bogoliubov vacuum state for the quasipar- 
ticles at hi. As the magnetic field is being turned off to 
zero, the state of the system \ip(t)} is getting excited from 
its instantaneous ground state. However, in a similar way 
as in Rcf. [5], we can follow the time-dependent Bogoli- 
ubov method and assume that the excited state \i>(i)) is 
a Bogoliubov vacuum for a set of time-dependent quasi- 
particle annihilation operators 

lm{t) = U* nm {t)c n + V nm {t)c n . (23) 

This Ansatz is a solution of Schrodinger equation when 
the Bogoliubov modes u nm (t) and v nm (t) solve time- 
dependent Bogoliubov-de Gennes equations 

. du n _ 

2h(t)u n - j n (U n+1 + V n +l) + jn-1 (v n -l - "n-l) , 
. dv n _ 
1 dt ~ 

-2h(t)u n + jn {v n +i + U n+ i) + jn-1 {V n -1 ~ U n -i) . (24) 

with the initial condition that at hi each mode is a pos- 
itive frequency eigenmode of the stationary BdG equa- 
tions (19). The evolution stops at h = when the final 
average density of kinks is 

d = a? E \ (! - IV'W) • (25) 

n 

Using the Jordan- Wigner transformation (14) followed 
by the Bogoliubov transformation (23), together with the 
assumption that 7(0)1^(0)) = 0, the density of kinks can 
be transformed into 

d = \~ >7 E R ° + "n+l,m(0)] < m (0) . 

ran 

(26) 
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The time-dependent BdG equations (24) were solved 
for a range of almost five decades of tq with the quench 
(22) starting at hi — 10 h c . Final densities of kinks were 
collected in Figure 2. To test the equation (9), valid 
for slow transitions, the data for the slowest quenches 
were fitted by (1{tq) — ^^ TQ ) w * tn ^ — 2. Here 
e{a,TQ) is a solution of Eq.(9) and a, (3 are two fitting 
parameters. Both a — 3.02 and (3 = 0.15 turned out 
to be 0(1) and the data are in reasonably good agree- 
ment with the best fit. Although Eq.(9) is confirmed by 
numerical simulations, even the slowest quenches with 
tq = 16384, obtained with substantial numerical effort, 
are not slow enough to clearly show the asymptotic be- 
havior in Eq.(ll) for tq — > oo. This is not unexpected 
because this asymptote is achieved when not only tq > 1 
but also In tq ^> 1. 

The same figure shows two power-law fits d <~ t q W ■ 
Relatively fast left-most quenches were fitted with the 
exponent w — 0.48 which is close to the 1/2 in the pure 
Ising model. These quenches, whose e is relatively large, 
cannot feel the randomness of J n 's and, as predicted, 
their final defect density scales like in the pure Ising 
model [5-7]. For comparison, the two right-most data 
points were also fitted by a power law d ~ Tq W , but this 
time the exponent is a mere w = 0.13. The exponent 
w decays with increasing tq as expected from the quasi- 
logarithmic solution e(a,TQ) of Eq.(9). 




Quench time 

FIG. 2. Density of kinks d as a function of transition time 
tq on a lattice of N = 512 sites. Circles are final kink densi- 
ties averaged over 4 realizations of J n . Error bars set the size 
of the circles. The solid line is the best fit <2(tq) = prgj[^yj 
where e(a,TQ) is a solution of Eq.(9). The dotted line is the 
best power law fit d ~ Tq W to the left-most 3 data points and 
the dashed line is the best power law fit to the right-most 2 
data points. The exponents are w = 0.48 and w — 0.13 re- 
spectively. The exponent w = 0.48 for the fastest transitions 
is consistent with the | in the pure Ising model. 



density of defects after a transition depends only loga- 
rithmically on the transition rate. Crudely speaking, the 
density is more or less the same no matter how slow is 
the transition. This result is in sharp contrast to the 
standard power law scaling predicted for pure systems. 
This strongly non-adiabatic behavior can be attributed 
to the anomalously small energy gap near the critical 
point which is characteristic for a random system with 
an infinite disorder fixed point. 
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Conclusion. — This paper is the first test of the 
Kibble-Zurek mechanism in a disordered system. The 
test is passed nicely but with an anomalous result that 



4 



